Draft version January 10, 201 1 

Preprint typeset using I^TgX style emulateapj v. 08/22/09 



A HIGH MASS DUSTY DISK CANDIDATE: THE CASE OF IRAS I8I51-1208 

C. Fallscheer' and H. Beuther 
Max Planck Institute for Astronomy, Konigstuhl 17, 691 17 Heidelberg, Germany 

J. Sauter^ and S. Wolf 

University of Kiel, Institute of Theoretical Physics and Astrophysics, Leibnizstrasse 15, 24098 Kiel, Germany 

AND 
Q. Zhang 

Harvard-Smithsonian Center for Astrophysics, 60 Gai'den St., Cambridge, MA 02138, USA 
Draft version January 10. 2011 

ABSTRACT 

Many questions remain regarding the properties of disks around massive prototstars. Here we present the 
observations of a high mass protostellar object including an elongated dust continuum structure perpendicular to 
the outflow. Submillimeter Array 230 GHz line and continuum observations of the high mass protostellar object 
IRAS 18151-1208 along with single dish IRAM 30m observations afford us high spatial resolution (0.8") as 
well as recovery of the extended emission that gets filtered out by the interferometer The observations of '-CO 
confirm the outflow direction to be in the southeast-northwest direction, and the 1.3 mm continuum exhibits 
an elongation in the direction perpendicular to the outflow. We model the physical parameters of the elongated 
structure by simultaneously fitting the observed spectral energy distribution (SED) and the brightness profile 
along the major axis using the 3D Radiative Transfer code MC3D. Assuming a density profile similar to that of 
a low mass disk, we can also reproduce the observations of this high mass protostellar object. This is achieved 
by using the same density distribution and flaring parameters as were used in the low mass case, and scaling up 
the size parameters that successfully modeled the circumstellar disk of several T Tauri stars. We also calculate 
that a region within the inner 30 AU of such a high mass disk is stable under the Toomre criterion. While we 
do not rule out other scenarios, we show here that the observations in the high mass regime are consistent with 
a scaled up version of a low mass disk. Implications on high mass star formation are discussed. 
Subject headings: stars: formation — stars: individual (IRAS 18151-1208) — stars: early type 



1. INTRODUCTION 

In contrast to low mass T Tauri stars and even intermediate 
mass Herbig stars of which numerous examples exist where 
circumstellar disks have been detected (see review by Wat- 
son et al. 2007), the role of disks in the case of massive star 
formation is not yet well defined. In the low mass case, it is 
clear that accretion disks and molecular outflows are integral 
parts of the formation process. While the occurrence of mas- 
sive molecular outflows is well established in cases of massive 
star formation, evidence for the analogous disk component 
has only tentatively been identified (see review by Cesaroni 
et al. 2007). 

In a concerted effort to shed some light on accretion disks in 
high mass star formation, we observed a sample of seven can- 
didates with the Submillimeter Array (Zhang et al. in prep). 
The high mass protostellar object IRAS 18151-1208 is one of 
the sources included in this study. This source was also in- 
cluded in the single dish survey of IRAS sources conducted at 
flie IRAM 30 m telescope by Sridharan et al. (2002); Beuther 
et al. (2002a), and has also been observed with the IRAM 
30 m and Mopra single dish telescopes by Marseille et al. 
(2008). IRAS 18151-1208 is at a distance of 3.0 kpc and has 
a luminosity of approximately 10"* L© (Sridharan et al. 2002). 

Previous single dish observations (Beuther et al. 2002b) 
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suggested the presence and general morphology of outflow 
in the region. However, data quality was not good enough to 
observe specific details. Narrow-band H2 {A = 2.122 fim) ob- 
servations (Davis et al. 2004) of IRAS 18151-1208 show evi- 
dence of two collimated molecular outflows originating from 
two separate sources. In the case of low mass star forma- 
tion, it is clear that accretion disks play a central role in the 
existence of outflows. However, in the high mass regime, ob- 
servations of scaled up versions of these low mass accretion 
disks are rare (Cesaroni et al. 2007). Our goal is to show that a 
scaled-up version of a disk modeling scheme that successfully 
reproduces the observations of disks around T Tauri stars can 
be applied in the case of a massive protostellar object, namely 
IRAS 18151-1208. 

The T Tauri star IRAS 04302 H-2247, nicknamed the "But- 
terfly Star" because of its butterfly-shaped morphology, has 
a wefl-defined circumstellar disk. Wolf et al. (2003, 2008) 
model the circumstellar environment of this source using a 
Monte Carlo 3D (MC3D) Radiative Transfer code. Similarly, 
the disk in the Bok globule CB 26 (Sauter et al. 2009) has also 
been observed and modeled. CB 26 is also a T Tauri star des- 
tined to become a low mass star with an edge-on circumstellar 
disk. Additionally, a handful of other edge-on disks around T 
Tauri stars have been modeled using a similar technique [e.g. 
HK Tau (Stapelfeldt et al. 1998), HV Tau (Stapelfeldt et al. 
2003), and IM Lupi (Pinte et al. 2008)]. On the theoretical 
side, Whitney et al. (2003) use a similar Monte Carlo radia- 
tive transfer approach to test various envelope/disk geome- 
tries for Class I sources. While theoretical descriptions such 
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Fig. 1. — The spectra from the SMA. Upper: the lower side band spectrum 
at the center of the 1.3 mm dust continuum peak of the combined compact 
and very extended configuration data. Middle: the same as above except for 
the upper side band. Lower: the upper side band at a position approximately 
7" northeast of the continuum peak only including the compact configuration 
data. 



as Whitney et al. (2003) and observational studies (Sicilia- 
Aguilar et al. 2006, for example) of low mass T Tauri disks 
are relatively well developed, such a complete picture is miss- 
ing in the high mass regime. 

2. OBSERVATIONS 

2.1. Submillimeter Array 

IRAS 18151-1208 was observed at 1.3 mm with the Sub- 
millimeter Array (SMA)^ in the compact and very extended 
configurations on 2007 July 8 and 2007 May 17 respectively. 
Observations in the compact array were made under stable 
and excellent weather conditions, with a t(230 GHz) of ap- 
proximately 0.07. The very extended configuration obser- 
vations were made under observing conditions with a t(230 
GHz) of 0.12. The combination of these two configurations 
provided baseline lengths varying between 14 and 520 meters 
corresponding to 10 k/l and 400 k/i at 1.3 mm. The phase ref- 
erence used was RA(J2000) = 18''17'"57.1' and Dec(J2000) 
— -12°07'22" with the emission peak approximately 17" east 
of these coordinates. We adopted a rest velocity Vlsr = 32.8 
km s ' (Sridharan et al. 2002). 

^ The Submillimeter Array is a joint project between the Smithsonian As- 
trophysical Observatory and the Academia Sinica Institute of Astronomy and 
Astrophysics, and is funded by the Smithsonian Institution and the Academia 
Sinica. 



For both data sets, the quasar 3C273 was used as the band- 
pass calibrator, and in between each 10 minute source obser- 
vation, 5 minute observations of the quasar 1743-038 were 
made for phase and amplitude calibration. The quasar 1911- 
201 was used as an additional phase and amplitude calibra- 
tor in the very extended data set. Fluxes for the calibra- 
tion sources were obtained from the Submillimeter Calibrator 
List"^. Measured fluxes after the calibration were consistent 
with the tabulated values to within 20%. A summary of this 
information can be found in Table 1 . 

For the mapping, a miriad robust parameter of was used 
for the continuum while a more natural weighting scheme 
(larger values of the robust parameter) was used for the line 
data. After combining the SMA data sets, the synthesized 
beam size of the continuum is 0.8"x0.7". A summary of 
this and other observation parameters are given in Table 1. 
Spectral resolution of 0.55 km s"' was obtained over a 4 GHz 
bandwidth divided into a lower and upper sideband separated 
by 10 GHz (see Fig. 1) with central frequencies of 219.6 and 
230.5 GHz respectively. The lines observed by the SMA are 
listed in Table 2. These are discussed in more detail in Section 
3.2. 

2.2. Pico Veleta 30 m Telescope 

The shortest baseline in our SMA datasets is 14 m which 
corresponds to a spatial scale of 19". The SMA is not sensi- 
tive to structures larger than this, so supplementary short spac- 
ing information was obtained in the on-the-fly mode with the 
HERA receiver on the IRAM 30 m telescope near Granada, 
Spain. These observations were made on 2006 November 12 
under very good observing conditions with t(230) of 0.07. 
The receiver was tuned to 230.5 GHz centered on the '^CO(2- 
1) line. A spectral resolution of 0.4 km s ' was obtained with 
a channel spacing of 0.3 MHz. The spectra were first pro- 
cessed with CLASS90 of the GILDAS package, then further 
analyzed with GREG. While the single dish data can in prin- 
ciple be combined with the interferometer data, their quality 
is not good enough to do so in this case, so we analyze the 
data sets independently. 

3. OBSERVATIONAL RESULTS 
3.1. Millimeter Continuum Emission 

Combining the very extended and compact configuration 
interferometer data, we map the 1.3 mm dust continuum at a 
spatial resolution of ~2500 AU (Fig. 2). Aside from the pri- 
mary peak, we detect a secondary source 16" to the southwest 
of the primary. Zooming into the primary peak (Fig. 3) dis- 
plays the elongation of the dust continuum in greater detail. 
This elongation is perpendicular to the bipolar molecular out- 
flow shown in Fig. 4, and may be associated with a disk com- 
ponent in a roughly edge-on orientation. The position angles 
of the relevant components are listed in Table 3. To simplify 
the analysis of the modeling, we have rotated the zoomed in 
image by 61° such that the elongation lies along the horizon- 
tal axis. The brightness profile along this axis is asymmetric, 
with the difference in height of the two peaks on the order of 
Icr (see Fig. 5). 

We measure an integrated flux of 0.4 Jy in the region con- 
tained within the 10,000 AU x 10,000 AU box shown in Fig. 3 
and a peak flux of 0.042 Jy. Following the methods of Hilde- 
brand (1983) and Beuther et al. (2002b, 2005) we adopt a 

http://smal.sma.hawaii.edu/callisl/callist.html 
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TABLE 1 

SMA OBSERVATION PARAMETERS. ENTRIES INCLUDE THE DATE OF OBSERVATION, THE ARRAY CONFIGURATION, THE CALIBRATORS, THE SYNTHESIZED BEAM OBTAINED AFTER 

INVERTING AND CLEANING THE DATA, AND THE IcT CONTINUUM NOISE LEVEL. 



Date 


Freq 
[GHz] 


Conliguration 


Calibrators 
Bandpass Phase & Amplitude 


beam 
["] 


rms„„, 
[mjy/bm] 


2007 May 17 
2007 Jul 08 


230 
230 


very extended 
compact 


3C273 1743-038 & 1911-201 
3C273 1743-038 








230 


combined 
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Fig. 2. — The 1.3 mm dust continuum. The size of the beam is indicated in the lower left corner. Contours start at 3cr increasing in steps of la where o" is 2 
mJy/beam. 
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Fig. 3. — A zoomed in view of the primary peak of the 1.3 mm dust contin- 
uum. Here we have rotated the image by 61° in order to align the elongation 
with the horizontal axis. The size of the beam is indicated in the lower left 
corner. The integrated flux in this area is 0.4 Jy. One arcsecond coiTesponds 
to 3000 AU. 



TABLE 2 
Observed lines 



Transition 


Rest Frequency 






[GHz] 


[K] 


C'»0 2 ^ 1 


219.560 


16 


SO 65 54 


219.949 


35 


'^CO 2 ^ 1 


220.399 


16 


CHiOH 8_i - 


7o E 229.759 


89 


C02^ 1 


230.538 


17 



TABLE 3 

Position angles (measured East of North) of the components associated 



WITH THE SOURCE. 


Component 


PA 




H2 jet 


128 


SMA '-CO outflow 


130 


30 m '^CO outflow 


108 


continuum elongation 


29 



dust opacity index ySop of 2.35 (Weingartner & Draine 2001) 
which corresponds to the dust model used in our modeling 
described in Section 4. Using a reference wavelength of 250 
yum, this opacity index corresponds to an opacity k of 0.19 
cm-g"' at 1.3 mm. Assuming a dust temperature of 30 K, 
and an optically thin system, our measured flux at 1.3 mm 
corresponds to a mass of 220 M© and a beam averaged col- 
umn density of 1.1x10^^ cm The mass and column density 
would of course be lower if we use different values for y6op, 
but to maintain consistency in the comparison with T Tauri 
stars, we use j6op of 2.35 for the modeling. Using the Inter- 
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stellar Medium value for j6op of 2.0 (/c=0.35 cm-g ') (Hilde- 
brand 1983), the mass and column density then would be 120 
Mo and 6.4 xlO^^ cm^^. With a value for ySop of 1.4 (/^=0.93 
cm^g ') (Ossenkopf & Henning 1994), the mass and column 
density become 45 Mq and 2.4 xlO^'^ cm"^. These masses 
are, of course, affected by missing flux in the interferometer 
data. Comparing our interferometric 1.3 mm flux to the single 
dish peak flux at this wavelength within the 11" beam of the 
MAMBO bolometer obtained by Beuther et al. (2002a), we 
estimate that less than 60% of the flux within this 11" beam 
is filtered out. 

3.2. Line Data 

Toward the continuum peak, there are only four spectral 
lines present in the 4 GHz bandwidth of the SMA data (see 
Fig. 1). Based solely on chemical evolution, it is difficult to 
determine whether this source has not yet reached the hot core 
phase, or whether it has already passed through this evolu- 
tionary phase. Three of the four lines present toward the con- 
tinuum peak, '2CO(2-l), '^€0(2-1), and SO(6-5), are not 
high density tracers, and C'^0(2-l) appears to be affected by 
the outflow (see Sect. 5.3). Hence, with the data at hand we 
are unable to infer the kinematics of the elongated structure 
detected in the continuum and discussed below. Methanol, 
a molecule linked to rotation in Fallscheer et al. (2009), is 
present a few (~7) arcseconds to the northeast of the primary 
peak but is extremely weak directly at the location of the pri- 
mary source. 

3.2.1. Outflow 

Previous single dish observations by Beuther et al. (2002b) 
were taken during non-ideal weather conditions, so the quality 
of our current data set is significantly improved. Nevertheless, 
the '^CO (2-1) emission is consistent with the outflow sug- 
gested by this older data set. The SMA data in Fig. 4 demon- 
strate that our observations follow the morphology of the H2 
jets seen by Davis et al. (2004). In both sets of observations, 
the two outflows observed are highly collimated. The two 
outflows we detect are roughly perpendicular to one another 
and appear to emanate from two separate sources detected in 
our 1.3 mm dust continuum. Outflows associated with low 
mass star formation depend on the presence of a disk. While 
massive outflows are likely also gas entrained by jets driven 
by magnetohydrodynamic acceleration from a disk, a descrip- 
tion of other methods is given in Arce et al. (2007). 

Several instances of anomalous behavior are apparent in 
Fig. 4. First, in the southeast-northwest outflow associated 
with the primary continuum source, there is a blueshifted 
component on the northwest side which is otherwise asso- 
ciated with redshifted emission (approximately located at 
AR.A.=4", ADec=7"). This might be an indication that the 
outflow has a small inclination angle relative to the plane of 
the sky such that we observe both receding and approaching 
sections of the outflow. Next, it is interesting to note that both 
sides of the outflow associated with the secondary continuum 
source are dominated by blueshifted emission. This region 
is dominated by redshifted emission in the single dish data 
(Fig. 4) indicating that it is large scale emission and is filtered 
out in the interferometric data. 

As discussed in Sect. 2.2, we do not combine the single dish 
data with the interferometer data. Instead, we use the interfer- 
ometer data to accurately determine the outflow morphology 
as described above, but rely on the single dish data to mea- 



sure fluxes. We use the single dish fluxes to calculate the out- 
flow mass, dynamical age and outflow energetics following 
the approach of Cabrit & Bertout (1990). Using an average 
spatial extent of 30" for each '^CO outflow lobe, and maxi- 
mum blueshifted and redshifted velocities of 10 and 17 kms"' 
respectively, we derive the properties listed in Table 4. Our 
derived values match those of Beuther et al. (2002b) reason- 
ably well, and are likely better estimates due to poor weather 
conditions at the time of the earlier observations. These val- 
ues confirm the original claim that the region will eventually 
form a high-mass star. 

3.3. Spectral Energy Distribution 

For the purpose of modeling, we gathered data for the spec- 
tral energy distribution (SED) of IRAS 18151-1208. The 
wavelengths and associated fluxes that we include in our SED 
(Fig. 6) are tabulated in Table 5. The 1.3 mm data point for 
the SED is determined by measuring the interferometric flux 
of IRAS 18151-1208 within the central 10,000 AU by 10,000 
AU region. We note that the flux measured with the inter- 
ferometer is a factor of four smaller than the corresponding 
integrated single dish flux reported by Marseille et al. (2008). 
The discrepancy between a factor 4 here and a factor of about 
2 discussed at the end of Sect. 3.1 comes from the fact that 
Marseille et al. (2008) integrate over a much larger area than 
just the 11" beam of the MAMBO observations. Since we 
would also like to use the 450 and 800 jjm fluxes reported in 
Marseille et al. (2008), we also divide these fluxes by a factor 
4 since they were measured over the same spatial scales as 
the 1.2 mm data. We justify this approach by the fact that the 
SED approximately follows a power law in this wavelength 
regime. 

For the SED, the 24.8 /im point comes from Campbell 
et al. (2008) and the 12.8 yum and 18.7 yum data points come 
from VISIR data (ESO proposal 075.C-0454, PI Fufler). The 
VISIR data were taken on 5 August 2006 with the VLT. The 
observations were made with the Q2 and [Neii] filters which 
have central wavelengths of 18.72 yum and 12.81 fim and 
widths of 0.88 and 0.21 /^m respectively. Integration time 
with the Q2 filter was 6 minutes and 5.5 minutes with the 
[Nen] filter. The standard star HD 1699 16 was used for flux 
calibration (Cohen et al. 1999). 

The nominal peak positions vary by ~1" between the mil- 
limeter and mid-infrared data. However, in the mid-infrared, a 
small field of view does not allow high positional accuracy so 
we attribute the emission to the same source. Alternatively the 
off'set may be due to a real difference in spatial position which 
would lead to a lowering of the SED at these wavelengths. We 
note this possibility, but proceed assuming that the detections 
at different wavelengths are spatially coincident. 

The data of Campbell et al. (2008) exhibit a silicate absorp- 
tion dip at 10 //m. This dip is likely caused by absorption 
of the envelope. Since the millimeter interferometer data fil- 
ter out the envelope, we are primarily interested in the small 
scale disk-like structure, so we do not try to reproduce this 
feature. 

4. MODELING 

After the 3D radiative transfer modeling of the circumstel- 
lar environment in the low mass regime was proven success- 
ful (Wolf et al. 2003; Sauter et al. 2009), we wanted to test 
whether a similar method could work in the high mass regime. 
We hypothesized that a scaled up version of a T Tauri disk 
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Fig. 4.— The outflow in ^^CO(2-\). Integrated maps of IRAS 18151-1208 over the entire velocity range of the line. Blueshifted emission (thick solid contours) 
integrated over velocities of 23-29 km s"' and redshifted emission (thick dashed contours) integrated over velocities of 40-50 km s"' . Left: SMA '^00(2-1) and 
continuum (grayscale with thin solid contours) overlaid on the H2 map of Davis et al. (2004). Contours for CO stall at 3(T and increase in steps of 3(T where cr is 
0. 15 and 0. 16 Jy bm"' -kms"' for the blue- and redshifted emission respectively. Right: IRAM 30m '-CO(2-l). The 1.3 mm dust continuum peak is indicated by 
the triangle, and the outflow direction and extent as determined by the SMA data are indicated by arrows. The region plotted in the SMA figure at left is indicated 
by the dotted lines. Contours start at 3o" and increase in steps of 3(T where cr is 11.4 Jy bm"' kms"' for the blueshifted component and 10.4 Jy bm"' kms"' for 
the redshifted component. 

TABLE 4 

Outflow properties derived from CO observations. Total outflow mass M|, momentum p, energy E, size, outflow dynamical age t, outflow rate Mom, 

MECHANICAL FORCE Fm AND MECHANICAL LUMINOSITY L,,, ARE GIVEN. INPUTS ARE DISCUSSED IN THE TEXT. 



M, [Mo] 


p [Mo km/s] E [erg] 


size [pc] t [yr] 


Mout [Mo/yr] 


F„, [Mo km/s/yr] 


Lm [L©] 


12 


190 3.1 X 10"" 


0.42 32000 


3.8 X K)-* 


6.0 X K)--' 


8.0 



TABLE 5 
Sources of SED data points. 



A 


iiux'' 


+/- 


aperture 


reference 


(41m) 


(Jy) 


(Jy) 


(") 




12.8 


26.3 


2 


5.0 


H. Linz, priv. comm. 


18.7 


41.0 


2 


5.0 


H. Linz, priv. comm. 


24.8 


101 


30 


1.9 


Campbell et al. (2008) 


450 


12.4 


3.7 




Marseille et al. (2008) 


800 


1.9 


0.6 




Marseille et al. (2008) 


1300 


0.4 


0.1 


3.3 


these data 



"The 450 and 800 pm fluxes are a factor of 4 smaller than the single dish 
fluxes pubUshed in Marseille et al. (2008) in order to make a more direct com- 
parison with the interferometer data. Further details are described in Sect. 3.3. 

may be able to reproduce the observations of what may pos- 
sibly be a disk around a massive protostellar object. 

While we are awai-e that the disk in IRAS 18151-1208 is 
still embedded in a large-scale envelope, we prefer to model 
only the disk structure. The interferometric map clearly shows 
an elongated structure, and while it may still have some con- 
tribution from a flattened envelope, applying a disk model to 
this structure appears to be a reasonable approach. Further- 
more, while we have missing flux on scales of the SMA pri- 
mary beam (55") on the order of 75%, it is only on the order 
of 60% when going to scales of the 11" beam of the IRAM 
30m observations (Sect. 3.1 and 3.3). Going to scales of the 
SMA synthesized beam (0.8"x0.7"), the missing flux is even 
less severe. Assuming that the single-dish peak flux within 



the 11" beam of 980 mJy (Beuther priv. comm.) is smoothly 
distributed, the single-dish flux per SMA beam is then 980 
mJy/[ll"xll"/(0.8"x0.7")] ~ 3 mJy which is only 1.5 times 
the rms of the SMA continuum image (see Tab. 1). Further- 
more, we can estimate the extinction at mid-infrared wave- 
lengths due to the envelope. Based on an extinction Ay of 72 
mag (Campbell et al. 2008), assuming Ry = 5.0 we estimate 
the extinction at 12 /^m to be approximately 2 mag (Mathis 
1990). Therefore, although a component is missing, model- 
ing only the disk appears to be a reasonable approximation for 
comparison purposes with low-mass disk models. 

4.1. The Monte Carlo 3D Radiative Transfer Code 

Using the self-consistent 3D Monte Carlo Radiative Trans- 
fer code, MC3D, developed by Wolf et al. (1999); Wolf 
(2003), we model the elongation seen in the dust continuum of 
IRAS 1 8 1 5 1 - 1 208 . In our setup, MC3D calculates the density 
distribution in a spherical region surrounding the central star 
divided into 91 equally spaced angular divisions and 50 log- 
arithmically spaced radial divisions such that the successive 
grid cell is ~1% longer than its inner neighbor This distri- 
bution is chosen to better resolve the density gradient in the 
more massive region closer to the central star. Using 100 log- 
arithmically distributed wavelengths between 0.05 and 20()0 
yum, MC3D computes the thermal structure of the disk, derives 
SEDs, and produces images at any of the 100 wavelengths. 
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For the purpose of this work, we use all 100 wavelengths for 
the temperature structure and the SED, and we compute an 
image at 1.3 mm. 

4.2. The Disk 

We use a density profile of a parameterized rotationally 
symmetric disk to produce fits to the observed SED and 
brightness profile along the elongation. To maintain a more 
even comparison, we assume the density profile used in the 
low mass cases (Wolf et al. 2003; Sauter et al. 2009) although 
we cannot exclude other density profiles with this approach. 
The density distribution we use is described by the following 
equation: 



Pdisk('"cyl) = PO ( — I exp 
v''cyl 



1 



(1) 



where z is the height above or below the midplane, R, is the 
radius of the central star and rcri is the radial distance from 
the central rotation (z-)axis. The normalization constant, po, 
is determined by the conditions placed on the total mass and 
volume of the system, and 



Ph = h 



100 



'"cyl 



lOOAU 



(2) 



where hun is the scale height of the disk 100 AU from the 
central star The shape of the modeled disk is described by the 
parameters a and /3 which are the radial density parameter and 
the disk flaring parameter respectively. If a and satisfy the 
relation a - 3(y6-|), then the density distribution corresponds 
to that of Shakura & Sunyaev (1973). The dependence of 
the surface density on the disk shape parameters can then be 
obtained by integrating Equation 1 with respect to z: £(rcyi) 
~ r^'yj where p-/3-a. 

4.3. The Dust Model 

Under the assumption that the gas is optically thin in the 
wavelength regime we are working in, it is the dust that is 
described by the density distribution in Eq. 1 . The dust grains 
are assumed to be spherical, and we adopt a power law size 
distribution following the MRN (Mathis et al. 1977) approach 
where the number of dust grains of a given radius a is given 
by 



«(fl) 



-3.5 



(3) 



Dust grain sizes a vary from 5 nm to 250 nm which are typical 
of the interstellar medium. 

We tested two different chemical compositions. The first 
follows the graphite and smoothed silicate mixture of Wein- 
gartner & Draine (2001) with relative abundances of 37.5% 
(graphite) and 62.5% (astronomical silicate). The other con- 
tains by mass 24.2% iron and iron-poor silicates, 5.6% troilite, 
25.6% refractory organics, 4.4% volatile organics, and 40.3% 
water (Pollack et al. 1994). 

4.4. The Parameter Space 

A motivating question behind this study is whether a low 
mass disk can be scaled up to describe a high mass disk. This 
would manifest itself in the form of having similar parameters 
describing the shape of the disk and scaled up parameters for 
the quantities such as disk mass and extent, accompanied by 
a more massive central star. In the scenario that scaling-up 



reasonably describes the observations of the high mass pro- 
tostellar object in question, we would expect the values for a 
and j6 to remain the same as in the low mass case. Although 
the value for hioo between the two cases is not directly compa- 
rable because 100 AU represents very different regions of the 
respective disks, we can scale hioo to a proportionally similar 
fraction of the disk's size and expect that to be proportionally 
larger than in the low mass case. In the ideal case, the disk 
mass should scale relative to the volume of the system which 
can be varied via the outer radius. 
The parameters that we include in the modeling are: 

• a, /3, and hioo: the density distribution exponent, the 
disk flaring exponent, and the scale height of the disk 
at 100 AU respectively. These quantities describe the 
dust density distribution via Equations 1 and 2. We test 
13 values of a between 1 and 4 and 1 1 values of /? be- 
tween 0.5 and 3. The relation a = 3(fi-\) is shown to 
hold in the case of low mass disks (Wolf et al. 2003; 
Sauter et al. 2009). We originally varied a and inde- 
pendently. However, after discovering that modifying 
a alone did not cause significant changes in the model- 
ing results, subsequent models followed the a - 3{/3-j) 
relation. We test 8 values of hioo between 1 and 30 AU. 

• 0: the inclination of the disk with respect to the plane 
of the sky. The observations indicate the presence of 
a flattened object, and this would only be visible at in- 
clination angles close to edge-on. However, we test 6 
inclination angles ranging from 45° to 90° (edge-on). 

• Rin: the inner radius. We initially chose small values 
for Rin of 10, 20, and 50 AU, but also tried significantly 
larger values such as 400, 1000, 2500, and 4000 AU 
before settling on 20 AU in our best fit model. 

• Rout: the radius at which the density goes to zero. The 
outer radius can roughly be constrained by the 1.3 mm 
continuum observations to a nominal value of 5,000 
AU, so we assume this value throughout. 

• Mdisk: the dust mass. Using the dust model of Weingart- 
ner & Draine (2001) we determine a rough estimate for 
the mass within the 5,000 AU radius to be 220 Mq (See 
details in Sect. 3.1). From this, we derive a dust mass of 
2.2 M0 based on a gas-to-dust ratio of 100. Taking this 
as our starting estimate, we subsequently determine the 
mass by ensuring that the SED fit goes through the 1.3 
mm data point. 

• Central star: due to the luminosity constraint of 
20,000 Lq described in the next section (Sect. 4.5), the 
most massive star we consider is a B0.5[T=26300 K, 
L=20000 La, R=7.0 Rq, M=15 Mq]. We also test a 
B1[T=25400 K, L=16000 Lq, R=6.4 Rq, M=13 Mq] 
and a B2 [T=22500 K, L=6000 Lq, R=5.6 Rq, M=10 
Mq] central star 

4.5. Observational Constraints 

For the SED, the fitting is constrained in a straightforward 
manner by the observed fluxes. We attempt to fit our SED over 
the entire wavelength regime for which we have data points 
available. 

While we use the SED to guide our efforts since the effects 
of modifying model parameters are more noticeable there (see 
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Figs. 7 and 8) than in the image, we aim to reproduce the elon- 
gated structure seen in the 1 .3 mm dust continuum perpendic- 
ular to the outflow orientation. Among other interpretations of 
what may be taking place in the region, possible explanations 
include a double source and a scaled up low mass disk. Here 
we focus on whether a scaled up low mass disk is consistent 
with these observations. Although we were able to reproduce 
a symmetric double peaked structure prior to smoothing with 
the PSF by making the inner radius larger, the convolution 
with the observed beam size smoothed out this effect in the 
image brightness profile. We then shifted our focus to fitting 
the shape of the elongated structure's outer edges. This is not 
a significant oversight because the difference in the height of 
the two peaks is only on the order of Icr. It is worth mention- 
ing that using an extremely large inner radius flattens out the 
peak of the brightness profile along the major axis (see Fig. 8), 
but because we do not know the origin of this slight asymme- 
try, we do not include anything in the modeling to account for 
the specific asymmetry nor the double peak and therefore do 
not give it any extra weight in our determination of the best fit 
model. 

For the model parameters, we were able to place rough con- 
straints on the central star, the outer radius, and the mass of the 
system before modeling began. Sridharan et al. (2002) place 
the luminosity of the IRAS 18151-1208 region at 20,000 Lq. 
We therefore use this as a limiting luminosity for the cen- 
tral star in our modeling. We use an outer radius of 5,000 
AU which is the radius at which the midplane density goes 
to zero and corresponds approximately to the 4cr contour of 
the elongation in the 1.3 mm dust continuum. Making the 
assumption that the region is optically thin, we can obtain a 
starting estimate for the mass of the system following the ap- 
proach described in Section 3.1. Our measured flux at 1.3 
mm corresponds to a mass of 220 Mq. We take this as a start- 
ing estimate, but since there are many sources of uncertainty 
in this calculation such as missing flux and unknown optical 
depth to name a few, we did not strictly limit the mass in our 
tests to this value. For the test cases with a gas mass of 220 
Mq, slight differences on the order of a few percent between 
the observed flux and the calculated flux at 1 .3 mm sometimes 
come about. We attribute this to the fact that despite using the 
same dust model, MC3D takes optical depth into considera- 
tion while the mass estimate does not. 

4.6. Comparing the Model with Observations 

The parameters of the best fit model are summarized in Ta- 
ble 6 and for comparison, the best fit model of the Butterfly 
Star (Wolf et al. 2003, 2008) and CB 26 (Sauter et al. 2009) 
are also included in this table. The corresponding SED and 
brightness profile of the best fit model are shown in Figs. 5 
and 6. All data points in the SED as well as the six points 
shown in Fig. 5 were used to determine the best fit model. 

We determine the best fit model based on considering both 
the SED and the brightness profile of the 1.3 mm image. For 
each test we calculate a reduced chi-squared value as de- 
scribed in the last paragraph of the next section (Sect. 5.1). 
The modeled image is convolved with the observed PSF and 
the brightness profile along the major axis is then obtained by 
summing the flux from the pixels in the columns perpendic- 
ular to the elongation. This is shown as the dashed line in 
Figs. 5 and 8. The pixels from columns aligned perpendicular 
to the elongation in the observed image are also summed and 
this is plotted as the solid line in Figs. 5 and 8. A comparison 
between the observed and modeled brightness curves is made 
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Fig. 6. — The SED of the best fit model. Due to differences in the observation 
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direct comparison with the 1.3 mm data point (see Sect. 3.3 for details.) 

by calculating the difference between the two curves at the six 
points indicated in Fig. 5. Weighting each point in the image 
and SED equally, we combine these differences with the dif- 
ference between the six observed and modeled SED points to 
arrive at the final value used to determine the best fit model. 
The numerical value of the weighted least squares calculation 
cannot be used on its own, but compared to the values from 
other models, we have a numerical value on which to base our 
decision of the best fit. Our best fit had a chi square value of 
2.4. While many models had chi square values less than 10, 
some were as high as 30. 

5. DISCUSSION 

5.1. Testing the Parameter Space 

As described in Sect. 4.4, a wide range of model parame- 
ters were tested before settling on the ones used in the best fit 
model. Figures 7 and 8 contain overlaid SEDs and brightness 
profiles from along the elongation that demonstrate the effects 
of varying the model parameters. The effects on the SEDs are 
conspicuous, but the differences in the brightness profiles are 
more subtle. We therefore quantify each fit by determining a 
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TABLE 6 

The parameters used m the best fit model. 
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Parameters include the radial density distribution parameter a, flaring parameter fS, scale height of the disk at jRout ho, dust mass Mjisi^, inclination angle 6 
(6=90° corresponds to an edge-on disk), the inner and outer radii of the density distribution Ri„ and Rom and the radius, temperature, and luminosity of the central 
star. 
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reduced chi-squared value for all tests performed. Although 
we do not have enough information about the envelope com- 
ponent to make quantitative statements about its effects, we 
discuss here how its omission affects our modeling results. 

As Figs. 7 and 8 demonstrate, some parameters can be con- 
strained better than others. We fix the outer radius. Rout based 
on the observations and although we still vary the values of 
the parameters for the disk's dust mass mjust and the lumi- 
nosity, temperature and radius of the central star, we do this 
mainly to see how varying these parameters affects the model. 
Ultimately, the values for these parameters were determined 
by the 1.3 mm observations rather than through the model- 
ing process. By not including an envelope component in the 
model, all of the millimeter dust emission is attributed to the 
disk. This means that the disk dust mass is likely overesti- 
mated. The outer radius as defined by the millimeter image 
may also be overestimated if the actual disk component is be- 
ing obscured by an envelope. The inner radius, on the other 
hand, is not well constrained by our models, and the addition 
of an envelope component would not affect the situation. 

The flaring exponent fi is constrained relatively well across 
the entire infrared and submillimeter wavelength regime as 
well as by the brightness profile, and we choose the corre- 
sponding value for the density distribution exponent a based 
on the relation known for low mass disks (Shakura & Sunyaev 
1973). Because some of the disk may be obscured by the en- 
velope component and especially the less dense outer flared 
regions, our value for the amount of flaring, characterized by 
j6, may be higher than if an envelope had been included in the 
modeling. Because of the fi - a relation, a would be similarly 
affected. However, the actual nature of the overestimation de- 
pends on the density distribution in the envelope, e.g. its radial 
power-law. 

The parameters for the disk scale height hioo, and the in- 
clination angle are predominantly constrained only by the 
SED in the mid-infrared wavelength regime. However, the in- 
clination angle is further limited by the millimeter image as 
it is seen as a flattened elongated object and hence is closer 
to edge-on opposed to face-on. As the envelope component 
plays more of a role at these shorter wavelengths, the val- 
ues derived for these parameters should be treated with some 
caution. In the presence of a thick spherically symmetric en- 
velope, the elongation of the disk would still be observed, 
providing support for the limitation on the inclination angle, 
but other envelope geometries may have different effects on 
the observed component. Using similar reasoning to the disk 
shape parameter the disk scale height would also be an 
overestimate because of the envelope component obscuring 
the outer part of the disk. 

Varying the dust composition as described in Section 4.3 
did not significantly affect our results, so we present here only 
the results of the Weingartner & Draine (2001) composition to 
maintain a more direct comparison with the Butterfly Star and 
CB 26. 

• The disk shape parameters fi (flaring parameter), and 
hioo (scale height at 100 AU) strongly affect the shorter 
wavelengths in the SED. As demonstrated in Fig. 7, 
both lower and higher values of fi and of hioo than 
used in the best fit model significantly decrease the 
flux in the shorter wavelength (micron) regime of the 
SED. Although the SED may be affected by a com- 
plex interplay between several effects, it is likely that 
the optical depth is primarily responsible for this be- 



havior Smaller values of beta (J3 < 1.0) mean that 
the disk is self-shadowed and thus less surface area is 
available for absorption and reemission of radiation. A 
more flared disk (larger value of P) obscures the central 
star more because of its larger surface area which inter- 
cepts more of the star's radiation. Both of these effects 
lead to a reduction of mid-infrared flux in the SED. For 
hioo, thicker disks (larger values) hide the central star 
while smaller values have the same amount of mass dis- 
tributed over a smaller volume leading to higher den- 
sity which increases the optical depth. As a result, the 
temperature gradient becomes steeper leaving a smaller 
area for reradiation which in turn means less flux in the 
micron regime. 

• The density distribution exponent a cannot be con- 
strained by the data themselves. Higher resolution data 
would be necessary to provide stronger constraints on 
a. After constraining [5, we use the direct relation for 
Shakura & Sunyaev (1973) disks between a and I3{a - 
lidS-^)} to determine a value of 2.4 for a. This approach 
satisfactorily reproduces our observations. 

• Inclination angle: the difference in the SEDs is small 
for inclination angles between 45° and 60°, but we use 
60° in our best fit model based on the observational data 
that indicate a more edge-on orientation. The SED de- 
viates further from the observations at inclination an- 
gles larger than 60°. At very high inclination angles 
(e.g. 80°), the upper layers of the disk block the central 
star from view and the SED drops significantly in the 
micron regime. 

• Disk mass: in the mass regime we tested, the only pa- 
rameter that significantly affects the millimeter wave- 
length side of the SED without affecting the micron side 
is the dust mass. Figure 7 makes it apparent that using 
the dust model of Weingartner & Draine (2001), a dust 
mass of 2.2 is necessary to fit the observed 1.3 mm 
flux. With a gas-to-dust ratio of 



100, 



(4) 



Dust 



this corresponds to a total disk mass of 220 M©. This is 
much larger than the main sequence mass of a Bl cen- 
tral star (13 Mq) indicating that this system cannot be 
Keplerian and should be quite unstable. To quantify this 
statement, we calculate the stability in the system un- 
der the assumption that it is a disk. Following Toomre 
(1964), the condition for stability is 



Q{r,z) 



nGUr) 



> 1 



(5) 



where Cj is the local sound speed, /Ck the angular fre- 
quency of the disk derived from Kepler's Law, G the 
gravitational constant and E the surface density. Fig- 
ure 9 indicates that regions only within the central-most 
30 AU of the system are stable under the Toomre crite- 
rion. These results are consistent with what Kratter & 
Matzner (2006); Vaidya et al. (2009) find in their mod- 
eling. As expected for a star-disk system that has a sig- 
nificant fraction of the mass in the disk, the modeled 
disk in this situation is unstable and is hkely fragment- 
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Fig. 8. — The effect on the 1.3 mm image's brightness profile of varying individual parameters. In each panel, the thick solid Hne is the cut through the midplane 
of the observed elongated structure, the dotted lines are +/- lo", the dashed line is the model and the thin solid line is an overlay of the best fit model (see Fig. 5). 
Here we have included only a select sample of the most deviant cases. The brightness profile was not affected as heavily as the SED and many of the cases show 
very minor deviations from the best fit model, so we do not include them here. In choosing the best fit model, we took both the SED and the brightness profile 
into account. See Sect. 4.6 for details. The parameters are similar to those of the best fit model except for the parameter(s) indicated in each panel. In panel {a): 
the disk flaring parameter /3=1.0 instead of 1.3. [b): the inclination angle 6=80° instead of 60° {c): the inner radius R|n=2500 AU instead of 20 AU and the dust 
mass Mdusi=1.2Mo instead of 2.2Mo. {d): the dust mass M(iust=1.2 Mq instead of 2.2 Mq. 

ing and possibly forming a multiple system as described 
by Krumholz et al. (2009). 

• Central star: because we do not have observational data 
for the wavelength regime dominated by the central 
star, we do not include these wavelengths in our SED 
modeling since the main contribution to the flux at the 
wavelengths for which we have data is from the warm 
dust in the disk. Nevertheless, increasing the temper- 
ature and luminosity of the central star raises the SED 
in the dust dominated wavelength regime and shifts the 
peak to shorter wavelengths as one would expect for 
warmer objects. 

• Inner radius: we do not expect our symmetric model 
to reproduce the asymmetric peak structure of the ob- 
served brightness profile along the elongation but even 
having a large inner hole (large Ri,,) does not produce a 
double peak structure as one might expect. At extreme 
values of Ri,, such as 2500 or 4000 AU, the bright- 
ness profile does appear double peaked, but after con- 
volution with the observed beam, the double peak that 
comes out of the modeling is smoothed out. This is 
not surprising given the ~2400 AU spatial resolution of 




200 



Fig. 9. — Contour plot of the Toomre Q stability parameter. The value of Q 
is labeled next to the respective contour. Values of Q greater than 1 .0 indicate 
stability under the Toomre criterion. 
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the data (see Table 1). Increasing the inner radius to 
these extreme values flattens out the brightness profile 
toward the center (see Fig. 8) but also makes the bright- 
ness profile wider which significantly worsens the fit at 
the edges. The SED is minimally affected by the in- 
ner radius except for very large inner radii since this 
removes hot dust close to the star that is responsible for 
mid-infrared flux. 

To quantify the quality of the fit for each set of parameters 
tested, we do a weighted linear least squares fit to the six data 
points in the SED as well as the six points in Fig. 5 distributed 
along the sides of the brightness profile. We avoid the central 
section of the brightness profile because our symmetric model 
cannot account for the asymmetric structure observed. We 
choose six approximately spatially independent points in the 
brightness profile to give equal weight to the SED and image 
fitting. This test helps significantly in narrowing down the 
best values for most parameters, although it does not do a 
good job discriminating between values of a. 

We have not done a complete parameter space study for the 
parameters listed above, so this least squares fit represents the 
best fit of the chosen parameter combinations. However, we 
have attempted to choose values for the parameters such that 
the progression of tests becomes more and more constraining 
for each parameter. Namely, starting with a coarse grid within 
the parameter space described in Sect. 4.4, a group of parame- 
ters was tested. By learning the effect on the SED and bright- 
ness profile of varying an individual parameter, subsequent 
tests of parameter sets were designed to target those parame- 
ters that would compensate for the most deviant aspects of a 
preceding model. This method narrowed down the parameter 
space and allowed finer sampling of the more realistic values. 

5.2. Comparison with T Tauri Stars 

Comparing IRAS 18151-1208 with CB 26 (Sauter et al. 
2009) and the Butterfly Star (IRAS 04302+2247, Wolf et al. 
2003, 2008), we notice that similar parameters describing the 
shape of the disk adequately reproduce the observed elonga- 
tion in IRAS 18151-1208. We did not choose these values 
right from the beginning, but only through quantifying the fits 
by means of calculating a weighted least squares value did we 
determine these parameters to be nearly identical. 

On absolute scales, the values of the disk scale height at 
100 AU do not coincide (5 AU for IRAS 18151-1208 versus 
10 AU for CB 26 and 15 AU for the Butterfly Star). How- 
ever, this discrepancy can be explained by the fact that in the 
T Tauri stars, 100 AU is proportionately much further out in 
the disk (ex. |Rout for the Butterfly Star) compared to the high 
mass protostellar object in which 100 AU corresponds to only 
one fiftieth of the outer radius (^Rout)- As defined, these val- 
ues cannot directly be compared. As a better comparison, we 
calculate the height of the massive disk at ^Rout- This turns 
out to be 200 AU, or 4% of the outer radius. In the CB 26 
and Butterfly Star cases, the height of the disk one third of 
the way out is 3% and 5% of the outer radius respectively. 
Within the uncertainties, these numbers are quite comparable 
and support our assumption that a geometrically thin disk is 
an appropriate model in the high mass regime. 

The parameters describing the dust mass and the central 
star differ significantly for the two types of young stars as ex- 
pected. The mass of the modeled IRAS 18151-1208 system 
is on the order of 1000 times greater than in the T Tauri cases. 
Although the calculated mass may be an overestimate if there 



is a contribution to the millimeter ffux from a dense envelope, 
we do not include an envelope in our models as that would 
increase the number of free parameters. 

Another large difference between the high mass protostellar 
object and T Tauri star systems is in the outer radius. In the 
scenario that massive star formation proceeds as a scaled up 
version of low mass star formation, it is expected that disks 
around massive stars should be significantly larger than the 
disks of their low mass counterparts. The observations in this 
case study are consistent with this picture, but definitive ob- 
servational evidence of large disks around massive stars is 
lacking (See the review by Cesaroni et al. 2007). The elon- 
gated structure that we observe in IRAS 18151-1208 is more 
than ten times larger than the several hundred AU extent of 
disks around T Tauri stars and although we do not want to 
overinterpret the scaling-up nature of our results, it is inter- 
esting to note that the disk's radius and the mass of the disk 
(which would increase by the cube of the scaling factor) can 
be described by a scaling factor of ~15 compared to the But- 
terfly Star. 

5.3. Rotation and Infall 

The scarcity of molecular lines in the observed spectrum 
has made it a difficult task to assess the region for indications 
of rotation or infall. A velocity gradient in C'^O is mildly 
present (see Fig. 10), but is neither aligned with the elonga- 
tion in the expected disk orientation nor with the outflow di- 
rection. This gradient has a velocity range of ~4 kms"' and 
extends across a region approximately 25,000 AU long. The 
SO emission follows a similar trend as C'**0, but the trend is 
even weaker These molecules are likely affected by both the 
outflow and the disk component. 

As for the other molecules detected in the SMA spectra (see 
Fig. 1), '"CO and '^CO are well established outflow tracers. 
The only other molecule detected in the region is methanol 
(CH3OH). However, there is no significant methanol detection 
at the position of the dust continuum peak, and the emission 
detected to the northeast does not exhibit any indication of 
rotation. As there are no further lines in our 4 GHz bandwidth 
spectrum from the SMA, we are limited on this front. Perhaps 
further observations at other wavelengths would yield a tracer 
of the kinematics in the region. 

5.4. Predictions for Future Observatories 

The MC3D code is capable of computing images at a wide 
range of wavelengths which can be used to speculate what a 
disk might look like with future observatories. In Fig. 1 1 we 
include predictions at 6, 12, and 25 jim as well as contour 
overlays of the 1.3 mm model. In this figure, we also include 
the beam sizes of observatories on the horizon, namely the 
Atacama Large Millimeter Array (ALMA), and the MlRl and 
METIS instruments planned for the James Webb Space Tele- 
scope (JWST) and the European Extremely Large Telescope 
(E-ELT). The flux is plotted on the same scale in each panel 
for the purpose of comparison. It should be noted that despite 
increases in sensitivity over today's instruments, the dynamic 
range of these future observatories will still be insufficient for 
the dynamic range of the 6 /^m image. 

6. CONCLUSIONS 

Observations of the high mass protostellar object IRAS 
18151-1208 reveal an elongation in the 1.3 mm dust contin- 
uum that has a position angle 29° East of North. These ob- 



IRAS 18151-1208 



13 



kms"^ 

31 32 33 34 35 36 
I I I I I I I I I I I I 



o 

O 
< 




32 34 
I ' ' ' I 




-5 - 



;t Moij|ent Map 




6000 AU 



31 32 33 34 35 36 
-| I I I I I I I I I I I I I I 



I I ' 1 TTT — I 1 1 — 

SO ist MomentnMap 



20 15 
AR.A. ["] 



25 



20 

ARA ["] 




20 15 
ARA ["] 



Fig. 10. — C'^'O, CH3OH, and SO velocity moment maps clipped at the 8, 4, and 60" levels respectively of the corresponding line's intensity map. For C"*0, 
iT=0.075 Jy, for CH3OH, (T=0.054, and for SO, o"=0.062. In each panel, the contours are those of the 1.3 mm dust continuum, the an'ows indicate the orientation 
of the '^CO outflow, and the eUipse at bottom left represents the size of the SMA beam. Although the size of each map is identical, note that the spatial position 
is offset in the middle panel compared to the other panels. 
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Fig. 1 1 . — Predictions at mid-infrared and millimeter wavelengths of what may be observable with future observatories. For the sake of example, resolution 
elements for several planned instruments (MIRl on JWST, METIS on the E-ELT, and ALMA) are shown in each panel that correspond to wavelengths observable 
with the respective instrument. The circles in the lower left corner con'espond to the beam sizes of the MIRI (white circles) and METIS (black circles) instruments 
and the white circle in the upper right corners represents a 0.1" ALMA beam. The flux scale is identical in all three panels. 



servations also confirm a well-defined orientation for a colli- 
mated massive outflow emanating from that dust continuum 
peak with a position angle 130° East of North. This elongated 
structure is perpendicular to the outflow orientation, and has 
an approximate diameter of 10,000 AU. While kinematic ev- 
idence of rotation is not available, the spatial signature of the 
dust continuum is indicative of a circumstellar disk. 

Under the hypothesis that massive star formation is a scaled 
up version of low mass star formation, we adapted the 3D 
radiative transfer Monte Carlo code MC3D to reproduce the 
observations. This was achieved by fitting the SED as well as 
the brightness profile along the major axis of the elongated 1.3 
mm dust continuum image. The best fit model is of a flared 
disk containing 220 of gas and dust with an inclination of 
60°. While this is not the only possible scenario to explain 
the observations, we conclude that this scaled up approach is 
viable. 

Compared to disks around T Tauri stars for which MC3D 
has successfully been applied, our best fit model has similar 
parameters describing the shape of the disk (density distribu- 
tion and amount of flaring) while larger values are necessary 



for the absolute parameters of the system such as mass, extent, 
and the central star. The disk produced by the modeling is un- 
stable aside from a small region close to the central star This 
is an expected outcome since the mass of the disk is dominant 
in the disk-star system. 

A further outcome of the modeling is that we can produce 
images at multiple wavelengths which can be helpful for visu- 
alizing what future observatories such as ALMA, JWST, and 
the ELTs will be capable of observing. Future work possible 
before these observatories become available includes looking 
for other tracers of the kinematics in the region. Pending the 
discovery of a molecule that traces rotation in the region, this 
source will be unique in that both the spatial density and flar- 
ing structure as well as the kinematic properties of the region 
will be known. 
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